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Summary 

Improvements that have been made to the COSMIC NASTRAN elements CIHEX.l and 
Q DM EMI are described. These elements are isoparametric representations of 
solid and membrane elastic behavior. Recent papers by the authors have 
shown the official COSMIC versions of these elements to be inferior to 
those available in the MacNeal-Schwendler Corporation (MSC) version of 
NASTRAN in that they are overly stiff for some loadings. Modifications 
have been made to these elements which reduce the order of integration for 
shear terms and, for the eight-mode solid element, add additional strain 
functions. The resulting element formulations give behavior similar to 
that of the MSC elements. The paper discusses the changes made in the 
element formulations and compares results of test problems with results 
from the official COSMIC elements and with the MSC elements. 


Introduction 


The isoparametric membrane quadrilateral element, QDMEM1, in COSMIC NASTRAN 
is a stand-alone element for use in modeling problems which exhibit plane 
stress behavior. It is a stand-alone element because there is no general 
plate element which currently uses the Q DM EMI for the membrane stiffness. 

In contrast, the MSC uses the ODMEM1 element for the membrane part of their 
QUAD4 general plate element. 

The results of a finite element idealization study using all of the 
available membrane elements in NASTRAN was reported in [1], Although the 
written version of the paper reported results only for elements available 
in MSC NASTRAN version 38, the version presented orally at the NASTRAN 
Colloquium showed results using both MSC-38 and COSMIC- 15.5. As presented 
at the colloquium, there was a marked difference in results for the QDMEM1 
elements from these two versions of NASTRAN. At the time, it was surmised 
that the discrepancy was due to a different manner in which the numerical 
integration was carried out in the two versions. In particular, it was 
shown that the COSMIC element exhibited overly stiff behavior for the 
problems investigated and that the MSC element was vastly superior. 

A similar study was conducted by the authors for solid elements and reported 
in [2]. This study was aimed at finding the best of the available solid 
elements to model thermal and gravity deformation effects on optical 
mirrors. Of particular importance was investigation of problems that might 
be encountered with elements that have aspect ratios in the range of 5 to 
10. Use of elements with this range of aspect ratio is necessary in 
modeling mirrors to avoid the need of extremely large models which would be 
required if element aspect ratios near unit were necessary. The results of 
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this study again clearly indicated the superiority of the MSC solid elements 
for modeling relatively thick plates for bending, as would occur for 
mirrors subjected to thermal gradients. For the eight-node solid isopara- 
metric elements, it was again shewn that there was a significant difference 
between the MSC element (HEXA-8) and the COSMIC element (CIHEX1). In 
addition, the COSMIC element showed extreme sensitivity to aspect ratio. 
Elements that had a thickness (in the plate thickness direction) smaller 
than its in-plane dimensions exhibited large errors. This was to be 
expected based on the above discussion of the membrane element deficiencies 
and the similarity of the formulation for the membrane and solid elements. 

The purpose of the effort reported herein, then, was to investigate modi- 
fications that could be made to the COSMIC Q DM EMI and CIHEX1 elements which 
would improve their accuracy for modeling bending type behavior. 

The next section discusses the cause of the overly stiff behavior of these 
elements. This has been investigated by others, [3] - [6], and is shown to 
be due to a parasitic shear that is introduced when these lower order 
isoparametric elements have bending modes of deformation. 

The following section describes the modifications that are required for the 
CIHEX1 element in order for it to behave as the MSC HEXA-8 element. The 
authors gratefully acknowledge the MSC for providing the mathematical 
description [7] of their modifications to the eight-node solid element and 
the QUAD4 general shell element. 

Following this is a brief description of the changes that are required for 
the QDMEM1 element, which consisted only of reducing the Gaussian integration 
order for the shear strain terms. 

Finally, the results of several problems are presented showing the improve- 
ment of the modified COSMIC elements in comparison to the officially 
installed element. 


Parasitic Shear in Bending 
in Lower Order Isoparametric Elements 

The quadrilateral membrane element, QDMEM1, has four grid points with 
stiffness for the 2 in-plane degrees of freedom at each grid point yielding 
a total of 8 degrees of freedom for the element. As discussed in [3], the 
8 degrees of freedom can be considered to be linear combinations of eight 
nodes of deformation, consisting of the three rigid body modes (two transla- 
tion and one rotation in its plane) plus the following deformation modes: 



(a) ( V) «) (d) 


x: location of Gauss Points for 2x2 integration 
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Modes (a) — (c) are the three constant strain modes while (d) and (e) are 
similar to bending modes. Generally, Gaussian quadriture is used to 
evaluate the stiffness matrix for isoparametric elements in which case the 
2x2 stiffness matrix for a pair of grid points is 

Cm s £** Sv 

The following section shows the details of this (for the 3-D element) with 
explanations of the terms. 

For now, it is sufficient to point out that the summation in the above 
equation is over the Gauss points (four in this case) and that C^i is the 
matrix relating displacements at grid point i to strains at Gauss point j. 
At these Gauss points, the terms in the C 9i matrix relating to shear 
strains are nonzero for the bending modes (d) and (e). Thus (d) and (e) 
modes will contribute shear strain energy in a situation where the element 
is used to model pure bending situations. In fact, as the element aspect 
ratio (i/d) increases, this parasitic shear becomes a dominant part of the 
strain energy and the element becomes excessively stiff for modeling 
bending. If, instead of evaluating the terms in which relate to shear 
strain at the Gauss points, terms were evaluated at the element center, it 
would be found that no shear strain energy would result in modes (d) and 
(e) since the shear strain is zero at the center. Then, since the shear is 
zero for modes (d) and (e) under this evaluation, modes (d) and (e) would 
indeed be pure bending; and the "effective" deformation in these modes 
would be: 



This is the motivation behind "reduced integration." One way to enforce 
the shear terms in to be evaluated at the center is to use a 

lxl Gaussian integration (f=l) in evaluating the above stiffness matrix 
equation. However, another way (which preserves the element volume) is to 
use the required Gaussian integration order needed to exactly evaluate the 
volume integrals (2x2 in this case) and evaluate the appropriate shear 
terms at the element center instead of at the Gauss point. It is this 
latter approach which is taken in modifying the COSMIC elements. 

The solid isoparametric element CIHEX1 has the same difficulty but in three 
dimensions and is modified with a similar reduced order integration for 
shear. 


Modified CIHEX1 Element 


For a 3-D solid isoparametric element, the displacements at any point in 
the interior of the element are expressed as (see Figure 9)): 
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( 1 ) 





where 



(?.) 


The summation on i is taken over all grid points of the element with A; 
being the vector of grid point i displacements. The are isoparametric 
interpolating functions in terms of the f coordinates which map the 
general hexahedron, in x,y,z coordinates into a rectangular parallelepiped 
inf ,T\ coordinates. For an eight-node hex element: 

Hi* j (3) 

with 5 1 the coordinates of grid point i 

The element strains are related to displacements by 


where 



(4) 
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is a vector of the element strains evaluated at the Gauss points 


and 
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(5) 


Note: ~ under a quantity indicates it is a matrix. 
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where the comma denotes partial differentiation with respect to x,y, or z 
and the subscript g on the matrix indicates that the terms in the matrix 
are to be evaluated for 5 at a Gauss point. This is the form of the 
strain-displacement relationship for the classical eight-node hex element 
and is what is employed in the COSMIC CIHEX1 element. Element stresses at 
Gauss points are related to strains through 

(5) 

where is the 6x6 constitutive matrix of material constants at a Gauss 
point. 

The element stiffness matrix is obtained using Gaussian quadriture from 

S» £« (7) 

where is a 3x3 partition of the element stiffness matrix relating 

forces at grid point i to displacements at grid point j. For exact inte- 
gration of the eight-node hex element, a 2x2x2 (g=8) Gaussian quadriture 
must be employed. Wg and Jg are Gauss integration weights and Jacobian. 

It has been observed that the stiffness matrix thus derived can exhibit 
overly stiff behavior when modeling bending situations due to the presence 
of parasitic shear in the element. One technique of overcoming- this diffi- 
culty is to employ a reduced integration (i.e., g less than 2x2x2). 

However, a general reduction of the integration order is not needed and can 
indeed lead to singular stiffness even when the element is restrained in a 
rigid body fashion. Instead, selective reduced integration is employed 
wherein only the terms in equation (7) which relate to shear strains have 
their integration order reduced. 

The HEXA eight-mode element in MSC NASTRAN employs reduced integration for 
shear. As explained by R. Harder [7], the MSC element uses for C-^the 
modified form: ~ 



32 



The terms in the upper half of equation (8) are identical to the form used 
in the classical isoparametric element. The terms in the lower portion of 
are different from their classical counterparts due to employing a 
reduced integration scheme. Harder proposes to use a Gaussian weighted 
average of the M • * , etc., to obtain , etc. 






T 




(9) 


with a similar definition for Ni jT . and Nj, . The summations in equation 
(9) are taken over some of the Gauss points: In particular, when evaluating 

terms in the fourth row of equation (8), the group summation is over all 
Gauss points in the plane in which / is a constant. For terms in the fifth 
row, the group is all points in a plane of constant J and, for the sixth 
row, a plane of constantly.. Harder shows that this Gaussian weighted 
averaging (thus reduced integration) scheme is necessary for the^element to 
maintain its capability to pass a constant strain patch test. 


In addition to employing this reduced integration for the shear terms, by 
averaging the related C^- M coefficients , Harder also employs additional 
"strain functions" in the MSC element to allow higher order polynomial 
variation of the direct strain terms. Strain functions are somewhat like 
bubble modes (see, for example, [8]), which have been used in some elements 
to also overcome their relatively stiff behavior in bending problems. 
Conceptually, the additional strain terms are included by modifying the 
basic strain-displacement relation of equation (4). Following the 
development in [7]: 

-•y r X Sja, + £}o (10) 

where is the classical c.. modified to represent the reduced order 

integration for shear. 


C<yois a 6xn matrix of strain coefficients and a. a vector of the amplitudes 
of the n strain functions added. Recognizing That the strain energy for a 
linear material element is 


Ui i Z V»sif» 

and minimizingUwith respect to the A 0 ampl itudes (keeping mind that is a 
function of through equations (6)*and (10), it is found that 


Using (6) and (10) 






in (11), the amplitudes are found as 


(ID 


& r K ' T K • & . 
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( 12 ) 
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where 



Finally, in combining (10) and (12), it is found that the strain displacement 
law in terms of only the physical degrees of freedom, is 


where 






£). 


(15) 

(16) 


With equation (15), the element stiffness matrix for the MSC eight-node hex 
element is generated as in equation (7) but with C,;, being used instead of 
C^;. It remains to select the strain functions to be added — i.e., the 
terms in the C, a matrix. In general, the C, 4 terms are added to fill a 
need in terms'of improvement in element accuracy for some particular 
application. Addition of these terms will almost certainly invalidate the 
interelement displacement continuity that exists with the classical 
element. However, this is not as significant as insuring that the modified 
element will still be capable of passing a constant strain patch test (see 
[9]) for a discussion of the patch test). Given the fact that the 
classical element does pass the constant strain patch test, Harder shows 
that the reduced integration technique will also pass the patch test and in 
order for the completely modified element to pass the patch test, it is 
required that: 

£».£**• < 17) 

For a constant patch stress, and keeping in mind that = 1.0 for the 
eight-node hex, equation (17) requires 


2 = o 


As pointed out in [7], this can be accomplished if the terms in 
the form _ . 0 

Thus, the MSC element uses for C^ 0 : 
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(18) 
are of 


(19) 
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Only terms are included in the first three rows as only the direct stresses 
are sought to be modified by the additional strain functions. It is not 
clear exactly why the particular form for was chosen, aside from the 
considerations in equation (18); however, it is the form identified by 
Harder in [7]. 

Equation (16) is the form of the strain-displacement law that was used in 
the modified CIHEX1 element reported herein, with equations (8), (13), 

(14), and (19) defining the various terms in- (16). This C.; matrix must 
be used in the development of the stiffness matrix as well as the thermal 
load vector and in stress data recovery. 


Modified QDMEM1 Element 


For the 2-D isoparametric element, the stiffness matrix has the same 
general form as shown in the previous section: 

?Uii 

However, for this element 
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Following the procedure of the previous section, reduced integration for 
shear terms is employed by modifying (22) to 
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( 24 ) 


As in the previous section, the M terms could be defined as Jacobian 
weighted averages of the terms at the Gauss points. However, for this 
element, it is found that the same result is obtained if the terns are 
evaluated at the center of the element ($=H=a). Thus 



MSC also uses strain functions for this element. From [7], for the 
membrane part of the QUAD4 general shell element 
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The complete formulation of C«; then follows as in the previous section and 
C«- would then be used in equation (20) instead of - For the modi- 
fications made to the COSMIC QDMEM1 element, however, the additional strain 
functions in were not implemented. Only the reduced integration for 

shear was implemented; and, therefore, C«. is used in equation (20) instead 


Code Changes to Implement Modified CIHEX1 Element 

The purpose of the authors' work was to prove the worth and feasibility of 
the element formulations. Optimizing the implementation of the code 
changes was not examined. Therefore, a complete discussion of the modi- 
fications will not be provided; however, a mention of the extent of the 
changes and the subroutines where they appear is appropriate. 

The alterations to the QDMEM1 were accomplished with a simple modification 
to the QDMM1D subroutine. At the point where' partitions of the stiffness 
are numerically integrated, terms which include the shear modulus were 
separated. This allowed the order of Gaussian integration to be selected 
separately for in-plane stress and shear stress terms. As discussed 
earlier, the need was for a linear shear stress formulation. The current 
implementation does not implement aniosotropic materials. 

The XIHEX subroutine, which calculates the mass and stiffness matrices for 
the CIHEX1 , as well as the CIHEX2 and CIHEX3 , received extensive modifica- 
tions. Since the IHEX subroutine recalculates part of the element's 
stiffness in order to form temperature loads, it required some similar 
changes. The CIHEX2 and CIHEX3 formulations were not considered, and 
alterations to the CIHEX1 stress recovery subroutine have not been developed 
at this time. 

The initial formulation of the CIHEX1 element matrices in COSMIC NASTRAfl 
was made in the basic coordinate system. Due to the manner in which the 
reduced integration for shear was implemented in the modified element, it 
was necessary to develop the CIHEX1 matrices in a local element coordinate 
system nearly aligned with the $ ,\,-J axes and then transform the matrices 
to the basic coordinate system. Additional subroutines needed to accomplish 
these tasks were written by the authors. Several options for selection of 
the "best" initial local system were looked at and are still being evalu- 
ated. For rectangular parallelepiped elements, the choice of a local 
system is trivial. For skewed elements, the? jX,'? directions do not form 
an orthogonal set of vectors in x,y,z; so the choice of the initial local 
coordinate system is not obvious. 
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Test Problems 


QDMEM1 Element 


The test problems are the same as those utilized in the prior mesh study of 
the membrane elements reported in [1] and involved a deep cantilevered beam 
type structure with unit depth and beam aspect ratio (length/depth) of two. 
Figure 1 shows the geometry, coordinate system, boundary conditions, and 
beam physical properties used in the study. Membrane elements formed the 
beam model. The mesh subdivision technique as well as the method used to 
indicate mesh size and element aspect ratio are demonstrated in 
Figure 2. 

The finite element model used work equivalent grid point forces for separate 

end moment and end shear loadings. This simulated the applied loads as 

well as the reactions at the cantilevered end (Figure 3). One should note, 
from Figure 1, that only kinematic constraints were imposed. Discussion of 
the theoretical solutions to these loading conditions can be found in [:]. 

In order to assess the effect of the reduced integration modification to 
the Q DM EMI element, the figures from [1] for the older membrane element 
mesh study were utilized. These curves, figure 4 through figure 8 herein, 
show the error in displacements and stresses at specific points on the beam 
as a function of mesh refinement or aspect ratio. The MSC element referred 
to in these figures is a CQUAD4 with only membrane properties specified on 
its PSHELL card; also, the 2,1 or 4,1 after COSMIC '84 refers to the number 

of Gaussian integration points for in-plane stress and shear stress terms, 

respectively. As indicated on the figures, the altered COSMIC element '84 
w/4,1 reduced integration produces the same answers as the MSC element; and 
these answers are an improvement over the old COSMIC element. 


CIHEX1 Element 

The premise for selecting a test problem in [2] was the fact that solid 
elements are used to model large optical mirrors of spaceborne telescopes, 
and these mirrors often have thicknesses of as much as 10 percent of their 
diameter. The test problem involved a cubic slab of equal dimension in the 
x-y plane and whose thickness varies between one-twentieth and one half of 
the x-y plane dimensions. Figure 10 shows the geometry, coordinate system, 
boundary conditions, and basic material information used in the study. The 
constraints are kinematic and the problem is symmetric about the x=0 plane. 
That is, the x displacement is zero along the x=0 plane. Using this 
constraint, only half the slab needed to be included in the finite element 
model. The mesh subdivision technique and method used to indicate element 
aspect ratio is shown in Figure 11. 

Originally, the test cases were chosen to measure the accuracy of various 
solid elements under temperature gradient and gravity loadings. Their 
value to this paper lies in that they provide separate bending and shear 
load cases. The linear temperature gradient produced a symmetric bending 
condition with a known theoretical answer. The gravity loading, hov/ever, 
was found to be non-converging and was used in this study only to show the 
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effect of the changes to CIHEX1 under a shear loading case. Further 
discussion of these load cases is provided in [2]. 

In order to assess effects of the changes to CIHEX1, three graphs were 
extracted from [2]. Data for the altered CIHEX1 was added to these graphs, 
which already included curves representing the original CIHEX1. 

Figures 13 and 14 present error in displacement, at a particular point, as 
a function of mesh size and aspect ratio, respectively. Actual displace- 
ments versus mesh size are shown in Figure 15 for the (non-convergent) 
gravity loading. Once again, the important factor to note in these graphs 
is that for each case the "improved" CIHEX1 provided the same answers as 
the comparable eight-node MSC element. For an aspect ratio of 10, not at 
all unreasonable when modeling large mirrors, the changes to the element 
totally eliminated a 48 percent error (in the temperature gradient case) 
when the old CIHEX1 element was used. 


Conclusions 


Modifications to the isoparametric membrane and solid elements, QDMEM1 and 
CIHEX1, have been implemented in the COSMIC NASTRAN code. 

The modified CIHEX1 element performs identically to the MSC HEXA eight-mode 
element. With the modifications, especially the reduced shear integration, 
it is anticipated that the new element will perform better when modeling 
thick plates when only few elements are used through the thickness. In 
addition, for pure bending, the element gives exact answers when only one 
element is used through the thickness. 

The modified QDMEM1 element has been shown to be superior to the original 
element when modeling bending situations. Neither element exhibits aspect 
ratio sensitivity in the modified form as it did in its original form. 
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NOTATION 


AR = PROBLEM ASPECT RATIO 

AR e = ELEMENT ASPECT RATIO 

ND = NUMBER OF ELEMENTS 
THROUGH DEPTH 


NL = NUMBER OF ELEMENTS 
ALONG LENGTH 


40 



FIG. 1 


BEAM GEOMETRY AND PROPERTIES 



£ = f .0508 M (2.0 IN) BASIC DEEP BEAM FOR MESH STUDY 
~ I VARIABLE FOR ASPECT RATION STUDY 

d= .0254 M (1.0 IN) 

t = .0254 M (1.0 IN) 

E= 1.9305 x 10 10 N/m 2 (28 x 10 6 LB/IN 2 ) 


v =0.3 


u = v = o AT x = y = o AND u = o AT x = o, y = —d/2 



FIG. 2 


FINITE ELEMENT MODEL MESH PATTERNS MEMBRANE ELEMENTS 


\» < H 



AR e = V d e 
AR = i/d 



TRMEM MESH PATTERN NO. 2 



| EXTRA NODE — QDMEM2 
I ELEMENT ONLY 
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TIP DEFLECTION ERROR, 


FIG. 4 

TIP DEFLECTION ERROR 
DEEP BEAM-END MOMENT LOADING 
(MESH SIZE STUDY) 



MESH SIZE, NL 




TIP DEFLECTION ERROR, 


FIG. 5 

TIP DEFLECTION ERROR 
DEEP BEAM-END SHEAR LOADING 
(MESH SIZE STUDY) 



MESH SIZE, NL 



DIRECT STRESS ERROR, 


FIG. 6 

DIRECT STRESS ERROR 

DEEP BEAM-END MOMENT AND END SHEAR LOADING 

(MESH SIZE STUDY) 



MESH SIZE, NL 




SHEAR STRESS ERROR, 


FIG. 7 

SHEAR STRESS ERROR 
DEEP BEAM-END SHEAR LOADING 
(MESH SIZE STUDY) 


33.3 



MESH SIZE, NL 



TIP DEFLECTION ERROR, % 


FIG. 8 

TIP DEFLECTION ERROR 
VARIABLE BEAM-END MOMENT LOADING 
(ASPECT RATIO STUDY) 




FIG. 9 

8 NODE HEX ELEMENT GEOMETRY 
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FIG. 10 

TEST PROBLEM-SOLID 
ELEMENTS 





♦ INDICATES DEGREE OF FREEDOM CONSTRAINED FOR 
KINEMATIC MOUNTS 

v = w = 0 AT Y = Z = 0, X = ±«/2 
u=w = 0ATX=Z = 0, Y =112 

A = (0.,0.,t) 

B = («/2,0.,t) 

MATERIAL INFORMATION (ALUMINUM) 

E = 6.89 x 10 10 N/M 2 (10 x 10 6 LB/IN 2 ) 
p =2.71 x 10 3 KG/M 3 (.098 LB/IN 3 ) 
a = 23.J x 10“ 6 /K 
u = .33 
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FIG. 11 


ELEMENT TYPES 



10 TETRAHEDRA (COSMIC) 

LINEAR ISOPARAMETRIC (COSMIC) 
LINEAR ISOPARAMETRIC (MSC) 


QUADRATIC ISOPARAMETRIC (COSMIC) 
QUADRATIC ISOPARAMETRIC (MSC) 


32 NODES 



CIHEX3 CUBIC ISOPARAMETRIC (COSMIC) 


CWEDGE 3 TETRAHEDRA (COSMIC) 
CPENTA LINEAR ISOPARAMETRIC (MSC) 


15 NODES CPENTA QUADRATIC ISOPARAMETRIC (MSC) 
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FIG. 12 


MESH GEOMETRY 



WEDGE PATTERN 
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Z DISPLACEMENT ERROR AT A,% 


FIG. 13 


Z DISPLACEMENT ERROR AT (A) 

THICK SLAB -LINEAR TEMPERATURE GRADIENT 

(MESH SIZE STUDY) 
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^CIHEXI (MODIFIED) 

ND = NL / 2 
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A = (0.,0.,t) 
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MESH SIZE, NL 




FIG. 14 

Z DISPLACEMENT ERROR AT (A) 

LINEAR TEMPERATURE GRADIENT 
(ASPECT RATIO STUDY) 
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Z DISPLACEMENT AT A, METERS 


FIG. 15 

Z DISPLACEMENT AT (A) 

THICK SLAB - GRAVITY 



MESH SIZE, NL 



